rm(list=ls());gc()
## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 534959 28.6 1193748 63.8 NA 669411 35.8
## Vcells 991680 7.6 8388608 64.0 16384 1851784 14.2
library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
## Warning: package 'ape' was built under R version 4.3.3
First, we need a relatedness matrix We created two phenotypic data set to discuss the role of envionmental and genetic (co)variances in G-matrices estimations
NB: in this tutorial we won’t ensure the convergence of the MCMC models or investigate autocorrelation in the posterior samples for simplicity
This tutorial has two parts
kin.inv=as(as.matrix(read.table('data/kinship.inverted.txt')),"dgCMatrix")
#dim(kin.inv) # 300
kin.inv[1:5,1:5]
## 5 x 5 sparse Matrix of class "dgCMatrix"
## ind1 ind2 ind3 ind4 ind5
## ind1 6.39710602 0.2796442 0.09781572 0.02587626 0.1746951
## ind2 0.27964425 6.8068225 -0.41487478 0.16668272 0.1517489
## ind3 0.09781572 -0.4148748 6.58310344 0.19790564 -0.1249038
## ind4 0.02587626 0.1666827 0.19790564 6.84654734 -0.3577153
## ind5 0.17469514 0.1517489 -0.12490376 -0.35771526 7.0852689
We will estimate the genetic variance for each trait separately and then in a single multivariate model
## Load the first data set
phen1 <- read.table(file="data/Pop1_phenotypes.txt",sep='\t',h=TRUE)
phen2 <- read.table(file="data/Pop2_phenotypes.txt",sep='\t',h=TRUE)
phen=phen1
#Rename the individuals to match the kinship matrix
phen$ind=row.names(kin.inv)
## Plot and investigate the phenotype distibution as you want
## Produce here your favorite plot
plot(phen$height,phen$weight,pch=16)
#cor(phen$height,phen$weight) # 0.6519626
We can see that the traits are highly correlated, but is it genetics or environment?
Run a first model on height
nb_trait=1
p.var = var(phen$height)
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height"))))
prior1.1 = list(G=list(G1=list(V=matrix(phen.var/2),nu=1)),R=list(V=matrix(phen.var/2),nu=1))
#model1.1 = MCMCglmm(height~1,random=~ind,ginv=list(ind = kin.inv),data=phen,prior=prior1.1,nitt= 53000, thin= 500, burnin= 3000)
#save('model1.1',file='RData/model1.1.RData')
load('RData/model1.1.RData')
# Convergence checks from Wednesday's workshop
heidel.diag(model1.1[["Sol"]]);heidel.diag(model1.1[["VCV"]])
##
## Stationarity start p-value
## test iteration
## (Intercept) passed 1 0.818
##
## Halfwidth Mean Halfwidth
## test
## (Intercept) failed -0.0134 0.0109
##
## Stationarity start p-value
## test iteration
## ind passed 1 0.691
## units passed 1 0.454
##
## Halfwidth Mean Halfwidth
## test
## ind passed 0.0942 0.00749
## units passed 0.9395 0.01522
autocorr.diag(model1.1[["Sol"]]);effectiveSize(model1.1[["Sol"]])
## (Intercept)
## Lag 0 1.00000000
## Lag 500 -0.03122956
## Lag 2500 -0.19051614
## Lag 5000 -0.04781202
## Lag 25000 -0.08813440
## (Intercept)
## 100
## Extract the genetic/residual variance estimates for the trait
plot(model1.1$VCV)
# Estimate genetic variance
apply(model1.1$VCV,2,median)
## ind units
## 0.08765121 0.93442541
# Credible interval for the trait
HPDinterval(model1.1$VCV)
## lower upper
## ind 0.04373865 0.1878074
## units 0.81512265 1.1066037
## attr(,"Probability")
## [1] 0.95
Run a second model on weight
nb_trait=1
p.var<-var(phen$weight)
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("weight"))))
prior1.2<-list(G=list(G1=list(V=matrix(phen.var/2),nu=1)),R=list(V=matrix(phen.var/2),nu=1))
#model1.2<-MCMCglmm(weight~1,random=~ind,ginv=list(ind = kin.inv),data=phen,prior=prior1.2,nitt= 53000, thin= 500, burnin= 3000)
#save('model1.2',file='RData/model1.2.RData')
load('RData/model1.2.RData')
plot(model1.2$VCV)
apply(model1.2$VCV,2,median)
## ind units
## 0.5801678 0.3939160
HPDinterval(model1.2$VCV)
## lower upper
## ind 0.4047552 0.8179593
## units 0.3196955 0.4844361
## attr(,"Probability")
## [1] 0.95
We see that two traits are highly phenotypically correlated, but models show one has very low heritability and one very high. We can run a multivariate model to estimate genetic (co)-variances between traits. So the environmental variance indeed is the reason behind phenotypic correlation.
Run a final model - multivariate
nb_trait=2
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height","weight"))))
prior1.3= list(G=list(G1=list(V=(phen.var/2),nu=2)),R=list(V=(phen.var/2),nu=2))
#model1.3<-MCMCglmm(cbind(height,weight)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,ginv=list(ind = kin.inv),data=phen,prior=prior1.3,
# family = rep("gaussian", nb_trait),nitt= 53000, thin= 500, burnin= 3000)
#save('model1.3',file='RData/model1.3.RData')
load('RData/model1.3.RData')
### Here we can extract a G-matrix of genetic (co)-variances
Gmat1 <- matrix(apply(model1.3$VCV[,1:4],2,median),nrow=2) # genetic co-variance matrix
# [,1] [,2]
#[1,] 0.10632641 0.08428513
#[2,] 0.08428513 0.45975364
matrix(apply(model1.3$VCV[,5:8],2,median),nrow=2) # phenotypic co-variance matrix (genetic + environment)
## [,1] [,2]
## [1,] 0.9371575 0.5748941
## [2,] 0.5748941 0.3970679
# [,1] [,2]
#[1,] 0.9371575 0.5748941
#[2,] 0.5748941 0.3970679
HPDinterval(model1.1$VCV)
## lower upper
## ind 0.04373865 0.1878074
## units 0.81512265 1.1066037
## attr(,"Probability")
## [1] 0.95
HPDinterval(model1.3$VCV)
## lower upper
## traitheight:traitheight.ind 0.06500519 0.1852113
## traitweight:traitheight.ind 0.01096517 0.1803158
## traitheight:traitweight.ind 0.01096517 0.1803158
## traitweight:traitweight.ind 0.33006258 0.5851530
## traitheight:traitheight.units 0.78990189 1.0881127
## traitweight:traitheight.units 0.49030694 0.6815341
## traitheight:traitweight.units 0.49030694 0.6815341
## traitweight:traitweight.units 0.33872217 0.4718267
## attr(,"Probability")
## [1] 0.95
Now we can compare the genetic/environmental variance estimated from both univariate and bivariate models and conclude on the source of phenotypic variance
Gmat shows genetic variance of trait1 (G11), trait2 (G12), and their co-varianc (G12) Similar for the environmental matrix, in the cor.test of 2 traits, we get R^2 = 0.66, this 0.66 is due to 0.57 environment and 0.08 only from genetics So the phenotypic correlation we got are actually environmental
Run similarly for Second population
phen = phen2
phen$ind=rownames(kin.inv)
### Explore the phenotypes as for population 1
plot(phen$height,phen$weight,pch=16)
#cor(phen$height,phen$weight) # 0.1711583
First model
nb_trait=1
p.var<-var(phen$height)
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height"))))
prior1.1<-list(G=list(G1=list(V=matrix(phen.var/2),nu=1)),R=list(V=matrix(phen.var/2),nu=1))
#model1.1_pop2<-MCMCglmm(height~1,random=~ind,ginv=list(ind = kin.inv),data=phen,prior=prior1.1,nitt= 53000, thin= 500, burnin= 3000)
#save('model1.1_pop2',file='RData/model1.1_pop2.RData')
load('RData/model1.1_pop2.RData')
plot(model1.1_pop2$VCV)
mean(model1.1_pop2$VCV)
## [1] 0.4934249
HPDinterval(model1.1_pop2$VCV)
## lower upper
## ind 0.3557060 0.8364216
## units 0.3129707 0.4819893
## attr(,"Probability")
## [1] 0.95
Second model
nb_trait=1
p.var<-var(phen$weight)
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("weight"))))
prior1.2<-list(G=list(G1=list(V=matrix(phen.var/2),nu=1)),R=list(V=matrix(phen.var/2),nu=1))
#model1.2_pop2<-MCMCglmm(weight~1,random=~ind,ginv=list(ind = kin.inv),data=phen,prior=prior1.2)
#save('model1.2_pop2',file='RData/model1.2_pop2.RData')
load('RData/model1.2_pop2.RData')
plot(model1.2_pop2$VCV)
mean(model1.2_pop2$VCV)
## [1] 1.189961
HPDinterval(model1.2_pop2$VCV)
## lower upper
## ind 0.5997141 1.722828
## units 1.0249987 1.542113
## attr(,"Probability")
## [1] 0.95
Bivariate model
nb_trait=2
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height","weight"))))
prior1.3<- list(G=list(G1=list(V=(phen.var/2),nu=2)),R=list(V=(phen.var/2),nu=2))
#model1.3_pop2<-MCMCglmm(cbind(height,weight)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,ginv=list(ind = kin.inv),data=phen,prior=prior1.3,
# family = rep("gaussian", nb_trait))
#save('model1.3_pop2',file='RData/model1.3_pop2.RData')
load('RData/model1.3_pop2.RData')
# Construct G/E matrices
#plot(model1.3_pop2$VCV)
Gmat2 <- matrix(apply(model1.3_pop2$VCV[,1:4],2,median),nrow=2)
# [,1] [,2]
#[1,] 0.3333886 -0.1420900
#[2,] -0.1420900 0.6488772
matrix(apply(model1.3_pop2$VCV[,5:8],2,median),nrow=2)
## [,1] [,2]
## [1,] 0.4434704 0.7345557
## [2,] 0.7345557 1.3751386
# [,1] [,2]
#[1,] 0.4434704 0.7345557
#[2,] 0.7345557 1.3751386
HPDinterval(model1.1_pop2$VCV)
## lower upper
## ind 0.3557060 0.8364216
## units 0.3129707 0.4819893
## attr(,"Probability")
## [1] 0.95
HPDinterval(model1.3_pop2$VCV)
## lower upper
## traitheight:traitheight.ind 0.2351921 0.46385776
## traitweight:traitheight.ind -0.2611148 0.02343029
## traitheight:traitweight.ind -0.2611148 0.02343029
## traitweight:traitweight.ind 0.3882450 0.96417018
## traitheight:traitheight.units 0.3665088 0.53121112
## traitweight:traitheight.units 0.6061634 0.87863519
## traitheight:traitweight.units 0.6061634 0.87863519
## traitweight:traitweight.units 1.1382431 1.68177154
## attr(,"Probability")
## [1] 0.95
We will create a null distribution forG matrices and discuss significance of the observed from the random
Generate a null distribution for the covariance The issue with MCMCglmm is that variance estimates are positive definite
## Randomize the individual labels
#phen.rdm <- phen
#all.post.modes <- NULL
#for(i in 1:100){
#phen.rdm$ind = sample(phen.rdm$ind)
#model1.3_pop2.rdm<-MCMCglmm(cbind(height,weight)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,ginv=list(ind = kin.inv),data=phen.rdm,prior=prior1.3,
# family = rep("gaussian", nb_trait))
#all.post.modes=rbind(all.post.modes,posterior.mode(model1.3_pop2.rdm$VCV))
#}
#save("all.post.modes",file='Rdm.model1.3_pop2.RData')
load("Rdata/Rdm.model1.3_pop2.RData")
# Get the credible interval from the randomized G matrices
apply(all.post.modes,2,function(x){HPDinterval(as.mcmc(x))})
## traitheight:traitheight.ind traitweight:traitheight.ind
## [1,] 0.08622752 -0.02542175
## [2,] 0.15828279 0.04040328
## traitheight:traitweight.ind traitweight:traitweight.ind
## [1,] -0.02542175 0.1750355
## [2,] 0.04040328 0.3195729
## traitheight:traitheight.units traitweight:traitheight.units
## [1,] 0.8652245 0.1739903
## [2,] 1.0038076 0.3098415
## traitheight:traitweight.units traitweight:traitweight.units
## [1,] 0.1739903 2.025782
## [2,] 0.3098415 2.279589
#traitheight:traitheight.ind traitweight:traitheight.ind
#[1,] 0.08622752 -0.02542175
#[2,] 0.15828279 0.04040328
#traitheight:traitweight.ind traitweight:traitweight.ind
#[1,] -0.02542175 0.1750355
#[2,] 0.04040328 0.3195729
#traitheight:traitheight.units traitweight:traitheight.units
#[1,] 0.8652245 0.1739903
#[2,] 1.0038076 0.3098415
#traitheight:traitweight.units traitweight:traitweight.units
#[1,] 0.1739903 2.025782
#[2,] 0.3098415 2.279589
Plot the intervals with the median of the estimated G-matrix and conclude
library(ggplot2)
# Observed G-matrix values
Gmat_observed <- c(0.3334, 0.6489, -0.1421)
# HPD intervals from randomized G-matrices
HPD_lower <- c(0.0862, 0.1750, -0.0254)
HPD_upper <- c(0.1583, 0.3196, 0.0404)
# Labels for each element of G-matrix
labels <- c("Var(height)", "Var(weight)", "Cov(height, weight)")
# Create dataframe for plotting
df <- data.frame(
Trait = labels,
Observed = Gmat_observed,
HPD_lower = HPD_lower,
HPD_upper = HPD_upper
)
# Plot the HPD intervals with observed values
# intervals show the intervals of HPD under random scenario (null distribution)
# the dots show the data point we observed
library(ggplot2)
ggplot(df, aes(x = Trait, y = Observed)) +
geom_point(size = 4, color = "salmon") + # Observed G-matrix values
geom_errorbar(aes(ymin = HPD_lower, ymax = HPD_upper), width = 0.2, color = "skyblue") +
theme_minimal() +
labs(y = "Genetic Variance/Covariance", x = "Matrix Element",
title = "Observed G-matrix vs. Null Distribution HPD Intervals") +
geom_hline(yintercept = 0, linetype = "dashed", color = "black")
We will then investigate how populations with different G matrices respond to selection NB: here we assume fix G-matrices that does not evolve only the mean population phenotype respond to selection
We can generate evolutionary path from the G-matrices The breeding value should align with the G matrix
## Plot the G-matrix as an ellipse
source("Rcode/99_Ellipse_function.R")
#plot(ellipse_coord(cor(phen1[,2:3])),xlim=c(-2,2),ylim=c(-2,2),type="l",col='skyblue') # P matrix
#lines(ellipse_coord(Gmat1)) # G matrix
#points(phen1$height/3,phen1$weight/3,pch=16,cex=.3)
library(ggplot2)
# Compute ellipse coordinates for P-matrix (phenotypic correlation) and G-matrix (genetic correlation)
P_ellipse <- as.data.frame(ellipse_coord(cor(phen1[, 2:3]))) # P-matrix
G_ellipse <- as.data.frame(ellipse_coord(Gmat1)) # G-matrix
# Rename columns for ggplot
colnames(P_ellipse) <- c("x", "y")
colnames(G_ellipse) <- c("x", "y")
# Scale the observed trait points for visualization
phen1_scaled <- data.frame(
x = phen1$height / 3,
y = phen1$weight / 3
)
# Plot the ellipses with ggplot
ggplot() +
geom_path(data = P_ellipse, aes(x = x, y = y), color = "skyblue", size = 1) + # P-matrix ellipse
geom_path(data = G_ellipse, aes(x = x, y = y), color = "chartreuse4", size = 1) + # G-matrix ellipse
geom_point(data = phen1_scaled, aes(x = x, y = y), color = "grey", alpha = 0.3, size = 0.5) + # Data points
theme_minimal() +
labs(title = "Comparison of G-matrix and P-matrix as Ellipses",
x = "Trait 1 (Height)",
y = "Trait 2 (Weight)") +
coord_fixed() # Ensure equal axis scaling
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Add more G matrix as needed
#create your own G matrix
Gmat3 = matrix(c(1,-0.3,-0.3,0.5),nrow=2)
G_ellipse2 = as.data.frame(ellipse_coord(Gmat3)) # G-matrix
colnames(G_ellipse2) = c("x", "y")
ggplot() +
geom_path(data = P_ellipse, aes(x = x, y = y), color = "skyblue", size = 1) + # P-matrix ellipse
geom_path(data = G_ellipse, aes(x = x, y = y), color = "chartreuse4", size = 1) + # G-matrix ellipse
geom_path(data = G_ellipse2, aes(x = x, y = y), color = "salmon", size = 1) + # second G-matrix ellipse
geom_point(data = phen1_scaled, aes(x = x, y = y), color = "grey", alpha = 0.3, size = 0.5) + # Data points
theme_minimal() +
labs(title = "Comparison of G-matrix and P-matrix as Ellipses",
x = "Trait 1 (Height)",
y = "Trait 2 (Weight)") +
coord_fixed() # Ensure equal axis scaling
We apply a vector of selection on 1 trait (directional selection)
beta1 = c(3,0) # the vector to "move" the matrix
#plot(ellipse_coord(Gmat1),xlim=c(-2,2),ylim=c(-2,2),type="l")
#new_mean_pop1 <- Gmat1%*%beta1
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta1
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta1
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
library(ggplot2)
# Generate ellipse coordinates for the G-matrix
G_ellipse = as.data.frame(ellipse_coord(Gmat1))
colnames(G_ellipse) <- c("x", "y")
# Compute mean shifts
new_mean_pop1 = Gmat1 %*% beta1 # move the matrix along the vector
mean_shifts = list()
for (i in 1:3) {
shifted_ellipse = as.data.frame(t(matrix(rep(new_mean_pop1, 1000), nrow=2)) + ellipse_coord(Gmat1))
colnames(shifted_ellipse) = c("x", "y")
mean_shifts[[i]] = shifted_ellipse
new_mean_pop1 = new_mean_pop1 + Gmat1 %*% beta1 # Update mean for next iteration
}
# Combine all shifts into a single dataframe
mean_shifts_df = do.call(rbind, Map(cbind, mean_shifts, shift_group = as.factor(1:3)))
# Plot using ggplot
ggplot() +
geom_path(data = G_ellipse, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_path(data = mean_shifts_df, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("salmon", "skyblue", "chartreuse4")) + # Different colors for each shift
theme_minimal() +
labs(title = "Genetic Covariance Ellipse with Mean Shifts",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed() # Keep aspect ratio fixed
The same with the second G-matrix
#plot(ellipse_coord(Gmat2),xlim=c(-2,2),ylim=c(-2,2),type="l")
#new_mean_pop2 <- Gmat2%*%beta1
#lines(t(matrix(rep(new_mean_pop2,1000),nrow=2))+ellipse_coord(Gmat2))
#new_mean_pop2 <- new_mean_pop2 + Gmat2%*%beta1
#lines(t(matrix(rep(new_mean_pop2,1000),nrow=2))+ellipse_coord(Gmat2))
#new_mean_pop2 <- new_mean_pop2 + Gmat2%*%beta1
#lines(t(matrix(rep(new_mean_pop2,1000),nrow=2))+ellipse_coord(Gmat2))
library(ggplot2)
# Generate ellipse coordinates for the G-matrix
G_ellipse2 = as.data.frame(ellipse_coord(Gmat2))
colnames(G_ellipse2) <- c("x", "y")
# Compute mean shifts
new_mean_pop2 = Gmat2 %*% beta1 # move the matrix along the vector
mean_shifts = list()
for (i in 1:3) {
shifted_ellipse = as.data.frame(t(matrix(rep(new_mean_pop2, 1000), nrow=2)) + ellipse_coord(Gmat2))
colnames(shifted_ellipse) = c("x", "y")
mean_shifts[[i]] = shifted_ellipse
new_mean_pop2 = new_mean_pop2 + Gmat2 %*% beta1 # Update mean for next iteration
}
# Combine all shifts into a single dataframe
mean_shifts_df = do.call(rbind, Map(cbind, mean_shifts, shift_group = as.factor(1:3)))
# Plot using ggplot
ggplot() +
geom_path(data = G_ellipse2, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_path(data = mean_shifts_df, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("salmon", "skyblue", "chartreuse4")) + # Different colors for each shift
theme_minimal() +
labs(title = "Genetic Covariance Ellipse with Mean Shifts",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed() # Keep aspect ratio fixed
We can generate any G-matrix and simulate the population response to selection
Gmat3 = matrix(c(1,-0.3,-0.3,0.5),nrow=2)
#plot(ellipse_coord(Gmat3),xlim=c(-5,5),ylim=c(-5,5),type="l")
#new_mean_pop3 <- Gmat3%*%beta1
#lines(t(matrix(rep(new_mean_pop3,1000),nrow=2))+ellipse_coord(Gmat3))
#new_mean_pop3 <- new_mean_pop3 + Gmat3%*%beta1
#lines(t(matrix(rep(new_mean_pop3,1000),nrow=2))+ellipse_coord(Gmat3))
# Generate ellipse coordinates for the G-matrix
G_ellipse3 = as.data.frame(ellipse_coord(Gmat3))
colnames(G_ellipse3) <- c("x", "y")
# Compute mean shifts
new_mean_pop3 = Gmat3 %*% beta1 # move the matrix along the vector
mean_shifts = list()
for (i in 1:3) {
shifted_ellipse = as.data.frame(t(matrix(rep(new_mean_pop3, 1000), nrow=2)) + ellipse_coord(Gmat3))
colnames(shifted_ellipse) = c("x", "y")
mean_shifts[[i]] = shifted_ellipse
new_mean_pop3 = new_mean_pop3 + Gmat3 %*% beta1 # Update mean for next iteration
}
# Combine all shifts into a single dataframe
mean_shifts_df = do.call(rbind, Map(cbind, mean_shifts, shift_group = as.factor(1:3)))
# Plot using ggplot
ggplot() +
geom_path(data = G_ellipse3, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_path(data = mean_shifts_df, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("salmon", "skyblue", "chartreuse4")) + # Different colors for each shift
theme_minimal() +
labs(title = "Genetic Covariance Ellipse with Mean Shifts",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed() # Keep aspect ratio fixed
Alternatively, we can define a phenotypic optimum and a gradient that linearly increase with phenotypic distance to the optimum
• Black ellipse → Original G-matrix covariance structure.
• Colored ellipses → Successive shifts in the mean phenotype.
#phen_optimum <- c(0,4)
#new_mean_pop1 <- c(0,0)
#plot(ellipse_coord(Gmat1),xlim=c(-.5,5),ylim=c(-.5,5),type="l")
#points(phen_optimum[1],phen_optimum[2],pch=16,col="red")
#beta_opt <- phen_optimum-new_mean_pop1
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta_opt
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#beta_opt <- phen_optimum-new_mean_pop1
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta_opt
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#beta_opt <- phen_optimum-new_mean_pop1
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta_opt
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
library(ggplot2)
# Define the phenotypic optimum and starting mean position
phen_optimum = c(0, 4)
new_mean_pop1 = c(0, 0)
# Compute the G-matrix ellipse
G_ellipse = as.data.frame(ellipse_coord(Gmat1))
colnames(G_ellipse) = c("x", "y")
# Store mean shift ellipses
mean_shifts = list()
# Compute the first shift
beta_opt = phen_optimum - new_mean_pop1
new_mean_pop1 = new_mean_pop1 + Gmat1 %*% beta_opt
# Compute and store 3 iterative mean shifts
for (i in 1:3) {
shifted_ellipse = as.data.frame(t(matrix(rep(new_mean_pop1, 1000), nrow=2)) + ellipse_coord(Gmat1))
colnames(shifted_ellipse) = c("x", "y")
mean_shifts[[i]] = shifted_ellipse
# Update new mean
beta_opt <- phen_optimum - new_mean_pop1
new_mean_pop1 <- new_mean_pop1 + Gmat1 %*% beta_opt
}
# Combine shifts into a dataframe
mean_shifts_df = do.call(rbind, Map(cbind, mean_shifts, shift_group = as.factor(1:3)))
# Convert phenotypic optimum into a dataframe for plotting
optimum_df = data.frame(x = phen_optimum[1], y = phen_optimum[2])
# Plot using ggplot2
ggplot() +
geom_path(data = G_ellipse, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_point(data = optimum_df, aes(x = x, y = y), color = "salmon", size = 4) + # Phenotypic optimum
geom_path(data = mean_shifts_df, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("skyblue", "chartreuse4", "darkgoldenrod1")) + # Different colors for each shift
theme_minimal() +
labs(title = "Evolutionary Mean Shifts Towards Phenotypic Optimum",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed(xlim = c(-0.5, 5), ylim = c(-0.5, 5)) # Set axis limits
Selection on trait 1 This plots a genetic covariance ellipse (G-matrix) and tracks mean shifts toward a defined phenotypic optimum. • The population mean starts at (0,0) and iteratively moves toward an optimum trait value at (0,4). • Each shift is calculated based on the difference between the current mean and the target optimum. • The shifts represent an adaptive trajectory where selection pushes the trait mean toward the optimum.
• Black ellipse → Original G-matrix covariance structure.
• Red point → Phenotypic optimum (target trait value).
• Colored ellipses → Successive mean phenotype shifts toward the optimum.
#phen_optimum <- c(4,0)
#new_mean_pop1 <- c(0,0)
#plot(ellipse_coord(Gmat1),xlim=c(-.5,5),ylim=c(-.5,5),type="l")
#points(phen_optimum[1],phen_optimum[2],pch=16,col="red")
#beta_opt <- phen_optimum-new_mean_pop1
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta_opt
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#beta_opt <- phen_optimum-new_mean_pop1
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta_opt
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#beta_opt <- phen_optimum-new_mean_pop1
#new_mean_pop1 <- new_mean_pop1 + Gmat1%*%beta_opt
#lines(t(matrix(rep(new_mean_pop1,1000),nrow=2))+ellipse_coord(Gmat1))
#... can be continued
library(ggplot2)
# Define the phenotypic optimum and initial mean position
phen_optimum = c(4, 0)
new_mean_pop1 = c(0, 0)
# Compute the G-matrix ellipse
G_ellipse = as.data.frame(ellipse_coord(Gmat1))
colnames(G_ellipse) = c("x", "y")
# Store mean shift ellipses
mean_shifts = list()
# Compute the first shift
beta_opt = phen_optimum - new_mean_pop1
new_mean_pop1 = new_mean_pop1 + Gmat1 %*% beta_opt
# Compute and store 3 iterative mean shifts
for (i in 1:3) {
shifted_ellipse <- as.data.frame(t(matrix(rep(new_mean_pop1, 1000), nrow=2)) + ellipse_coord(Gmat1))
colnames(shifted_ellipse) <- c("x", "y")
mean_shifts[[i]] <- shifted_ellipse
# Update new mean
beta_opt <- phen_optimum - new_mean_pop1
new_mean_pop1 <- new_mean_pop1 + Gmat1 %*% beta_opt
}
# Combine shifts into a dataframe
mean_shifts_df <- do.call(rbind, Map(cbind, mean_shifts, shift_group = as.factor(1:3)))
# Convert phenotypic optimum into a dataframe for plotting
optimum_df <- data.frame(x = phen_optimum[1], y = phen_optimum[2])
# Plot using ggplot2
ggplot() +
geom_path(data = G_ellipse, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_point(data = optimum_df, aes(x = x, y = y), color = "salmon", size = 4) + # Phenotypic optimum
geom_path(data = mean_shifts_df, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("skyblue", "chartreuse4", "darkgoldenrod1")) + # Different colors for each shift
theme_minimal() +
labs(title = "Evolutionary Mean Shifts Towards Phenotypic Optimum (Trait 1 Selection)",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed(xlim = c(-0.5, 5), ylim = c(-0.5, 5)) # Set axis limits
For pop2 and any pop similarly
#phen_optimum <- c(0,4)
#new_mean_pop2 <- c(0,0)
#plot(ellipse_coord(Gmat2),xlim=c(-.5,5),ylim=c(-.5,5),type="l")
#points(phen_optimum[1],phen_optimum[2],pch=16,col="red")
#beta_opt <- phen_optimum-new_mean_pop2
#new_mean_pop2 <- new_mean_pop2 + Gmat2%*%beta_opt
#lines(t(matrix(rep(new_mean_pop2,1000),nrow=2))+ellipse_coord(Gmat2))
#beta_opt <- phen_optimum-new_mean_pop2
#new_mean_pop2 <- new_mean_pop2 + Gmat2%*%beta_opt
#lines(t(matrix(rep(new_mean_pop2,1000),nrow=2))+ellipse_coord(Gmat2))
library(ggplot2)
# Function to compute mean shifts and return ellipse data
compute_mean_shifts <- function(G_matrix, phen_optimum, steps = 3) {
mean_shifts <- list()
new_mean <- c(0, 0) # Initial population mean
for (i in 1:steps) {
beta_opt <- phen_optimum - new_mean
new_mean <- new_mean + G_matrix %*% beta_opt
shifted_ellipse <- as.data.frame(t(matrix(rep(new_mean, 1000), nrow=2)) + ellipse_coord(G_matrix))
colnames(shifted_ellipse) <- c("x", "y")
mean_shifts[[i]] <- cbind(shifted_ellipse, shift_group = as.factor(i)) # Store shift step
}
return(do.call(rbind, mean_shifts))
}
# Define phenotypic optimum
phen_optimum <- c(0, 4)
# Compute ellipses for each population
G_ellipse2 <- as.data.frame(ellipse_coord(Gmat2))
colnames(G_ellipse2) <- c("x", "y")
mean_shifts2 <- compute_mean_shifts(Gmat2, phen_optimum)
G_ellipse3 <- as.data.frame(ellipse_coord(Gmat3))
colnames(G_ellipse3) <- c("x", "y")
mean_shifts3 <- compute_mean_shifts(Gmat3, phen_optimum)
# Convert phenotypic optimum to dataframe
optimum_df <- data.frame(x = phen_optimum[1], y = phen_optimum[2])
# Plot for Population 2
p2 = ggplot() +
geom_path(data = G_ellipse2, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_point(data = optimum_df, aes(x = x, y = y), color = "salmon", size = 4) + # Phenotypic optimum
geom_path(data = mean_shifts2, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("skyblue", "chartreuse3", "darkgoldenrod1")) + # Different colors for each shift
theme_minimal() +
labs(title = "Population 2: Evolutionary Mean Shifts",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed(xlim = c(-0.5, 5), ylim = c(-0.5, 5)) # Set axis limits
# Plot for Population 3
p3 = ggplot() +
geom_path(data = G_ellipse3, aes(x = x, y = y), color = "black", size = 1) + # G-matrix ellipse
geom_point(data = optimum_df, aes(x = x, y = y), color = "salmon", size = 4) + # Phenotypic optimum
geom_path(data = mean_shifts3, aes(x = x, y = y, group = shift_group, color = shift_group), size = 1) + # Shifted ellipses
scale_color_manual(values = c("skyblue", "chartreuse3", "darkgoldenrod1")) + # Different colors for each shift
theme_minimal() +
labs(title = "Population 3: Evolutionary Mean Shifts",
x = "Trait 1",
y = "Trait 2",
color = "Shift Step") +
coord_fixed(xlim = c(-2, 5), ylim = c(-2, 5)) # Set axis limits
p2
p3
Backsolving the breeding values to estimate SNP effect
We can use the genotype matrix to estimate SNP effects from breeding
values Caution: GBLUP are shrunk (meaning extreme values are
underestimated?) because of ridge regression (REML) or uninformative
priors, so SNP estimates are not accurate using this method
load('RData/Backtransformation_matrices.RData')
## We need to run model1.3 again and keep the posterior distribution of the
## random effects by adding pr=TRUE
#nb_trait=2
#phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height","weight"))))
#prior1.3<- list(G=list(G1=list(V=(phen.var/2),nu=2)),R=list(V=(phen.var/2),nu=2))
#model1.3_pop2_with_rand.post<-MCMCglmm(cbind(height,weight)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,ginv=list(ind = kin.inv),data=phen,prior=prior1.3,
# family = rep("gaussian", nb_trait),nitt= 53000, thin= 500, burnin= 3000,pr=TRUE)
#save('model1.3_pop2_with_rand.post',file='RData/model1.3_pop2_with_rand_post.RData')
load('RData/model1.3_pop2_with_rand_post.RData')
col_subset = colnames(model1.3_pop2_with_rand.post$Sol)[grepl(paste0("weight.ind."), colnames(model1.3_pop2_with_rand.post$Sol))]
breeding_values = data.frame(model1.3_pop2_with_rand.post$Sol[, col_subset])
new_colnames = gsub(".*ind\\.", "", colnames(breeding_values))#
colnames(breeding_values) = new_colnames
# Convert breeding values to matrix
breeding_values_matrice = as.matrix(breeding_values)
## data frame of chr position
# I'm just making up a random SNP matrix here
genotypes = read.csv("data/Simulated_Genotype_Data.csv", row.names=1)
#only need 200 SNPs
genotypes = genotypes[,1:200]
snps_positions = data.frame(chr='I',pos=1:ncol(genotypes))
# Calculate SNP effects
SNPs_effects.weight = M_transpose %*% MM_prime_inverse %*% t(breeding_values_matrice)
## And for height
col_subset <- colnames(model1.3_pop2_with_rand.post$Sol)[grepl(paste0("height.ind."), colnames(model1.3_pop2_with_rand.post$Sol))]
breeding_values <- data.frame(model1.3_pop2_with_rand.post$Sol[, col_subset])
new_colnames <- gsub(".*ind\\.", "", colnames(breeding_values))#
colnames(breeding_values) <- new_colnames
# Convert breeding values to matrix
breeding_values_matrice <- as.matrix(breeding_values)
SNPs_effects.height <- M_transpose %*% MM_prime_inverse %*% t(breeding_values_matrice)
snps_positions$weight <- apply(SNPs_effects.weight,1,median)
snps_positions$height <- apply(SNPs_effects.height,1,median)
### Load the SNP effect as generated
plot(snps_positions$height)
true_geno_effects=read.table("data/Pop2_geno_effects.txt",h=TRUE)
plot(snps_positions$weight,true_geno_effects$height)
summary(lm(snps_positions$height~true_geno_effects$weight))
##
## Call:
## lm(formula = snps_positions$height ~ true_geno_effects$weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.090753 -0.020568 0.002637 0.026268 0.082732
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.645e-05 2.593e-03 0.026 0.98
## true_geno_effects$weight -1.309e-01 2.715e-02 -4.820 2.86e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03631 on 198 degrees of freedom
## Multiple R-squared: 0.105, Adjusted R-squared: 0.1005
## F-statistic: 23.23 on 1 and 198 DF, p-value: 2.856e-06
summary(lm(snps_positions$weight~true_geno_effects$height))
##
## Call:
## lm(formula = snps_positions$weight ~ true_geno_effects$height)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.095351 -0.027590 0.003952 0.029384 0.098740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001077 0.003044 -0.354 0.724
## true_geno_effects$height -0.183195 0.028973 -6.323 1.67e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04294 on 198 degrees of freedom
## Multiple R-squared: 0.168, Adjusted R-squared: 0.1638
## F-statistic: 39.98 on 1 and 198 DF, p-value: 1.67e-09
We will here split the G/E matrices using within/between line variance instead of using a pedigree or a kinship matrix (see previous workshops)
Load the first data set
phen = read.table(file="data/Pop1_phenotypes_second_WS.txt",sep='\t',h=TRUE)
head(phen)
## ind height weight
## 1 1 -2.2445172 -0.66770736
## 2 2 -0.4380629 0.02265420
## 3 3 1.1791095 -0.12778455
## 4 4 -1.3845064 0.04073303
## 5 5 -0.7754505 -0.36484391
## 6 6 -0.8109581 -0.27744043
NB: some redundancy with the first workshop… but some students did not attend them Re-introduce rapidly the animal model / Introduce the Price equation
Run first an animal model with the two phenotypic traits Run a final model - multivariate
nb_trait=2
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height","weight"))))
prior1.3 = list(G=list(G1=list(V=(phen.var/2),nu=2)),R=list(V=(phen.var/2),nu=2))
#model1.3<-MCMCglmm(cbind(height,weight)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,data=phen,prior=prior1.3,
# family = rep("gaussian", nb_trait),nitt= 53000, thin= 500, burnin= 3000)
#save('model1.3',file='RData/model1.3_WS2.RData')
load('RData/model1.3_WS2.RData')
Gmat1 = matrix(apply(model1.3$VCV[,1:4],2,median),nrow=2)
Same thing for the population 2 that we will use later
phen2 = read.table(file="data/Pop2_phenotypes_second_WS.txt",sep='\t',h=TRUE)
head(phen2)
## ind height weight
## 1 1 -1.1285808 -0.08775198
## 2 2 0.4538317 -0.18213257
## 3 3 -0.1416177 -2.62585728
## 4 4 -0.1575402 1.69044269
## 5 5 -0.5649287 -0.37458336
## 6 6 -1.2021183 -0.55627078
nb_trait=2
phen.var = diag(nb_trait) * diag(var(subset(phen2, select = c("height","weight"))))
prior1.3 = list(G=list(G1=list(V=(phen.var/2),nu=2)),R=list(V=(phen.var/2),nu=2))
#model1.3.pop2<-MCMCglmm(cbind(height,weight)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,data=phen2,prior=prior1.3,
# family = rep("gaussian", nb_trait),nitt= 53000, thin= 500, burnin= 3000)
#save('model1.3.pop2',file='RData/model1.3_pop2_WS2.RData')
load('RData/model1.3_pop2_WS2.RData')
Gmat2 = matrix(apply(model1.3.pop2$VCV[,1:4],2,median),nrow=2)
Generate fitness estimates for the individuals - add some noise Initial scenario where only one trait is under selection
w = 4-phen$height
phen$w=w/mean(w)
## STNS - estimate a G matrix with fitness - to predict phenotypic evolution
nb_trait=3
phen.var = diag(nb_trait) * diag(var(subset(phen, select = c("height","weight","w"))))
prior1.4 = list(G=list(G1=list(V=(phen.var/3),nu=3)),R=list(V=(phen.var/3),nu=3))
#model1.4<-MCMCglmm(cbind(height,weight,w)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,data=phen,prior=prior1.4,
# family = rep("gaussian", nb_trait),nitt= 53000, thin= 500, burnin= 3000)
#save('model1.4',file='RData/model1.4.RData')
load('RData/model1.4.RData')
Gzw = matrix(apply(model1.4$VCV[,1:9],2,median),nrow=3)
print(Gzw)
## [,1] [,2] [,3]
## [1,] 0.034491080 0.04453024 -0.007685951
## [2,] 0.044530238 0.54472688 -0.011167551
## [3,] -0.007685951 -0.01116755 0.002155797
###############################################################################
### Here we need to build a random estimates: are the estimates larger than ##
### expected by chance ? ##
###############################################################################
phen.rand = phen
#post.Gzw <- NULL
#prior1.4<- list(G=list(G1=list(V=(phen.var/3),nu=3)),R=list(V=(phen.var/3),nu=3))
#for(i in 1:100){
# print(i)
#phen.rand$ind=sample(phen.rand$ind)
#model1.4.rand<-MCMCglmm(cbind(height,weight,w)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,data=phen.rand,prior=prior1.4,
# family = rep("gaussian", nb_trait),verbose=FALSE)
#post.Gzw <- rbind(post.Gzw,apply(model1.4.rand$VCV,2,median))
#}
#save('post.Gzw',file='RData/post.Gzw.Rand_model4.RData')
load("RData/post.Gzw.Rand_model4.RData")
HPDinterval(as.mcmc(post.Gzw)[,c(3,6,9)])
## lower upper
## traitw:traitheight.ind -0.008342513 -0.005449856
## traitw:traitweight.ind -0.004603998 -0.002157472
## traitw:traitw.ind 0.001572915 0.002295653
## attr(,"Probability")
## [1] 0.95
#lower upper
#traitw:traitheight.ind -0.008342513 -0.005449856
#traitw:traitweight.ind -0.004603998 -0.002157472
#traitw:traitw.ind 0.001572915 0.002295653
apply(model1.4$VCV[,c(3,6,9)],2,median)
## traitw:traitheight.ind traitw:traitweight.ind traitw:traitw.ind
## -0.007685951 -0.011167551 0.002155797
#traitw:traitheight.ind traitw:traitweight.ind traitw:traitw.ind
#-0.007685951 -0.011167551 0.002155797
Recapitulate the selection gradient Now that we have the G-matrix and the expected phenotypic change over 1 generation, we can deduce the gradient
beta_coefs = solve(Gmat1)%*%Gzw[1:2,3]
print(beta_coefs)
## [,1]
## [1,] -0.165201736
## [2,] -0.005679371
plot(phen$weight,phen$w)
summary(lm(height~w,data=phen))
##
## Call:
## lm(formula = height ~ w, data = phen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.080e-16 -2.130e-16 -5.800e-17 1.000e-16 1.777e-13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.000e+00 2.462e-16 1.624e+16 <2e-16 ***
## w -4.000e+00 2.389e-16 -1.674e+16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.271e-15 on 2998 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.803e+32 on 1 and 2998 DF, p-value: < 2.2e-16
summary(lm(weight~w,data=phen))
##
## Call:
## lm(formula = weight ~ w, data = phen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.58349 -0.54852 -0.00311 0.53547 2.52359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.57348 0.05765 44.64 <2e-16 ***
## w -2.57348 0.05593 -46.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7657 on 2998 degrees of freedom
## Multiple R-squared: 0.4139, Adjusted R-squared: 0.4137
## F-statistic: 2117 on 1 and 2998 DF, p-value: < 2.2e-16
And so we can build a credible interval for the beta coefficients
beta_distribution = NULL
for(i in 1:nrow(model1.4$VCV)){
beta_distribution = cbind(beta_distribution,solve(matrix(model1.4$VCV[i,1:9],nrow=3)[1:2,1:2])%*%Gzw[1:2,3])
}
HPDinterval(as.mcmc(t(beta_distribution)))
## lower upper
## var1 -0.33413444 -0.144381136
## var2 -0.01630024 0.009251241
## attr(,"Probability")
## [1] 0.95
Generate different evolutionary scenario -> generate the Gzw matrices and see how it matches the expectations ?
We have created different fitness distribution for the population 2, you need to analyze them, predict the phenotypic evolution and deduce the selection on the different phenotypes
phen2$w=read.table("data/fitness_pop2.txt",h=TRUE)$x
nb_trait=3
phen.var = diag(nb_trait) * diag(var(subset(phen2, select = c("height","weight","w"))))
#prior1.4_pop2 <- list(G=list(G1=list(V=(phen.var/3),nu=3)),R=list(V=(phen.var/3),nu=3))
#model1.4_pop2<-MCMCglmm(cbind(height,weight,w)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,data=phen,prior=prior1.4_pop2,
# family = rep("gaussian", nb_trait))
#save('model1.4_pop2',file='RData/model1.4_pop2.RData')
load('RData/model1.4_pop2.RData')
Gzw_pop2 = matrix(apply(model1.4_pop2$VCV[,1:9],2,median),nrow=3)
#phen2.rand = phen2
#post.Gzw_pop2 <- NULL
#prior1.4<- list(G=list(G1=list(V=(phen.var/3),nu=3)),R=list(V=(phen.var/3),nu=3))
#for(i in 1:100){
# print(i)
#phen2.rand$ind=sample(phen2.rand$ind)
#model1.4.rand2<-MCMCglmm(cbind(height,weight,w)~trait-1,random=~us(trait):ind,
# rcov=~us(trait):units,data=phen2.rand,prior=prior1.4,
# family = rep("gaussian", nb_trait),verbose=FALSE)
#post.Gzw_pop2 <- rbind(post.Gzw_pop2,apply(model1.4.rand2$VCV,2,median))
#}
#save('post.Gzw_pop2',file='RData/post.Gzw_pop2.Rand_model4.RData')
load('RData/post.Gzw_pop2.Rand_model4.RData')
min(angle_theta(e2[,1],EV_gamma[,2]),angle_theta(e2[,2],-EV_gamma[,2])) # 15°
## [1] 15.65939
min(angle_theta(eigen(cor(phen[,2:3]))$vectors[,1],EV_gamma[,1]),angle_theta(eigen(cor(phen[,2:3]))$vectors[,1],-EV_gamma[,1])) # 15°
## [1] 15.48188
source("Rcode/99_Ellipse_function.R")
#plot(phen2[,2],phen2[,3],col=adjustcolor(colors[color_index],alpha=.2),pch=16,cex=.6)
#lines(ellipse_coord(cor(phen[,2:3])),col='darkgreen',lwd=2)
#lines(ellipse_coord(Gmat2),col='orange',lwd=2)
# Convert color indexing to a valid color format with transparency
phen2$plot_color <- scales::alpha(colors[color_index], alpha = 0.2)
# Compute ellipse coordinates
P_ellipse <- as.data.frame(ellipse_coord(cor(phen[, 2:3]))) # P-matrix ellipse
G_ellipse <- as.data.frame(ellipse_coord(Gmat2)) # G-matrix ellipse
colnames(P_ellipse) <- c("x", "y")
colnames(G_ellipse) <- c("x", "y")
# Create ggplot
ggplot() +
geom_point(data = phen2, aes(x = height, y = weight, color = plot_color), size = 0.6, alpha = 0.2) + # Scatter plot
geom_path(data = P_ellipse, aes(x = x, y = y), color = "cyan3", size = 0.6) + # P-matrix ellipse
geom_path(data = G_ellipse, aes(x = x, y = y), color = "salmon", size = 0.6) + # G-matrix ellipse
theme_minimal() +
labs(title = "Phenotypic and Genetic Covariance Ellipses",
x = "Trait 1",
y = "Trait 2",
color = "Color Index") + scale_color_brewer(palette = "RdYlBu") + # Color scale
coord_fixed() + # Maintain aspect ratio
guides(color = "none") # Remove the color legend
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette RdYlBu is 11
## Returning the palette you asked for with that many colors
## Warning: Removed 131 rows containing missing values or values outside the scale range
## (`geom_point()`).